////////////////////////////////////////////////////////////////
// Sunil Suram (PhD Candidate)								  //
// Prof. Mark Bryden                                          //
// Virtual Engineering Research Group                         //
// Iowa State University, Ames, IA, 50010                     //
////////////////////////////////////////////////////////////////

#include "pod.h"

////////////////////////////////////////////////////////////////////////////////
pod::pod()
{
   //set default interpolation type to linear
   interpolationType = 1;
}

////////////////////////////////////////////////////////////////////////////////
pod::~pod( )
{

	gsl_matrix_free( coeffs );
	gsl_vector_free( vectorOfMeanData );
	gsl_matrix_free( meanRemovedData );

	if(mode=="construct")
	{
		gsl_matrix_free( inputData );
		gsl_matrix_free( v );
		gsl_vector_free( s );
		gsl_matrix_free( copyOfMeanRemovedData );
	}

	if(mode=="eval")
	{
		gsl_vector_free( PODapproximation );
	}
}
////////////////////////////////////////////////////////////////////////////////
double pod::Interpolation( const vector<double>& val, 
      const multi_array<double,2>& out )
{
	double ev = 0;
	BiCubicInterpolation* cubic = new BiCubicInterpolation();
	//cubic->SetInputPoints
	return ev;
}
////////////////////////////////////////////////////////////////////////////////
double pod::Interpolation( const vector<double>& x, 
      const vector<double>& out )
{
	double ev = 0;
	//inter->SetInputPoints
	return ev;
}

////////////////////////////////////////////////////////////////////////////////
double pod::Interpolation( const vector<double>& val, 
   const multi_array<double,3>& out )
{
   double ev = 0;
   if ( interpolationType == 2 )
   {
      cout<<"Performing linear interpolation of coefficients "<<endl;
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 3 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateCubicSplineAt( val, out );
      cout<<"eval val : "<<ev<<endl;
   }
   if ( interpolationType == 1 )
   {
      cout<<"Performing cubic spline interpolation of coefficients "<<endl;
      
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 3 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateLinearAt( val, out );
      cout<<"eval val : "<<ev<<endl;
      
   }
   
   return ev;

}
////////////////////////////////////////////////////////////////////////////////
double pod::dist( vector<double> v1, vector<double> v2 )
{
	//distance between 2 vectors, assumes the sizes of the vectors are the same and goes with
	//the vector v1
	double d = 0, r1=0, r2=0;
	for ( size_t i=0;i<v1.size();i++ )
	{
		r1 = r1 + v1.at(i)*v1.at(i);
		r2 = r2 + v2.at(i)*v2.at(i);
	}
	r1 = sqrt( r1 );
	r2 = sqrt( r2 );
	d = fabs( r1 - r2 );
	return d;	
}
////////////////////////////////////////////////////////////////////////////////
double pod::length( vector<double> v1 )
{
	double r1 = 0;
	for ( size_t i=0;i<v1.size();i++ )
		r1 = r1 + v1.at(i)*v1.at(i);
	return ( sqrt( r1 ) );
}
////////////////////////////////////////////////////////////////////////////////
double pod::Interpolation( const vector<double>& val, 
   const multi_array<double,4>& out )
{
   double ev = 0;
   if ( interpolationType == 2 )
   {
      cout<<"Performing linear interpolation of coefficients "<<endl;
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 4 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateCubicSplineAt( val, out );
      cout<<"eval val : "<<ev<<endl;
   }
   if ( interpolationType == 1 )
   {
      cout<<"Performing linear interpolation of coefficients "<<endl;
      
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 4 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateLinearAt( val, out );
      cout<<"eval val : "<<ev<<endl;
      
   }
   
   return ev;

}
double pod::Interpolation( const vector<double>& val, 
   const multi_array<double,5>& out )
{
   double ev = 0;
   if ( interpolationType == 2 )
   {
      cout<<"Performing linear interpolation of coefficients "<<endl;
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 5 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateCubicSplineAt( val, out );
      cout<<"eval val : "<<ev<<endl;
   }
   if ( interpolationType == 1 )
   {
      cout<<"Performing cubic spline interpolation of coefficients "<<endl;
      
      // generate interpolated coefficients
      MultiVariateInterpolation* mv = new MultiVariateInterpolation( );
      mv->SetNumberOfParameters( 5 );
      cout<<"set params :"<<endl;
      mv->SetInputPoints( inputParameterMap, out );
      cout<<"set inputs :"<<val.size()<<endl;
      ev = mv->EvaluateLinearAt( val, out );
      cout<<"eval val : "<<ev<<endl;
      
   }

   return ev;
}
////////////////////////////////////
//////PRIVATE MEMBER FUNCTIONS//////
////////////////////////////////////



////////////////////////////////////////////////////////////////////////////////
void pod::SubtractMeanOfData( )
{
	//this function computes the mean of the data

	//and subtracts the mean from the original data set
	vectorOfMeanData = gsl_vector_calloc( sizeOfDataSet );

	//compute mean now
	double tmpMean;
	FILE *foutMean;
	foutMean=fopen( ".\\Results\\mean.txt","wb" );

	for ( int j=0;j<sizeOfDataSet;j++ )
	{
		tmpMean = 0;
		for ( int i=0;i<numDataSets;i++ )
		{
			tmpMean = tmpMean + gsl_matrix_get( inputData, j, i );
		}
		tmpMean = tmpMean/(double)(numDataSets);
		gsl_vector_set( vectorOfMeanData, j, tmpMean );
	}
	gsl_vector_fprintf(foutMean,vectorOfMeanData, "%g");
	fclose(foutMean);

	//subtract mean of data from original data
	double tmp;
	for ( int i=0;i<numDataSets;i++ )
	{
		for ( int j=0;j<sizeOfDataSet;j++ )
		{
			tmp = gsl_matrix_get( inputData, j, i ) - 
               gsl_vector_get( vectorOfMeanData, j );
			gsl_matrix_set( meanRemovedData, j, i, tmp );
			gsl_matrix_set( copyOfMeanRemovedData, j, i, tmp );
		}
	}
}

////////////////////////////////////////////////////////////////////////////////
void pod::ComputeSVD( )
{
	/////////////////////////////////////////////////////////////////////////////
	//compute the singular value decomposition (SVD) of the mean removed data  //
	//A = USV'                                                                 //
	//A is the mean removed data matrix                                        //
	//U is the matrix of left singular vectors                                 //
	//S is the vector of singular values                                       //
	//V is the matrix of right singular vectors                                //
	//' is the transpose operator                                              //   
	/////////////////////////////////////////////////////////////////////////////

	//This procedure computes the <thin SVD> of the given matrix
	//Currently using this GSL based SVD computation as a starting point
	//might have to implement the full SVD algorithm from Golub/VanLoan book 
   //"Matrix Computations"

	if ( DBG_POD ) cout<<"Computing SVD ..."<<endl;
	//gsl_linalg_SV_decomp_jacobi( meanRemovedData, v, s );
	gsl_vector* work = gsl_vector_calloc( numDataSets );
	gsl_linalg_SV_decomp( meanRemovedData, v, s, work );

	//normalize singular values
	double tmp;
	for ( unsigned int i=0;i<s->size;i++ )
	{
		tmp = gsl_vector_get( s, i )/gsl_vector_get( s, 0 );
		if ( DBG_POD ) cout<<tmp<<endl;
	}

	//if ( DBG_POD > 1 )
	{
		FILE* foutU = fopen ( ".\\Results\\uMatrix.txt", "wb");
		gsl_matrix_fprintf( foutU, meanRemovedData, "%g" );
	}

	gsl_vector_free( work );

	//if ( DBG_POD > 1 )  //write matrices/vectors to files if you need them
	{

		FILE* foutV = fopen ( ".\\Results\\vMatrix.txt", "wb");
		gsl_matrix_fprintf( foutV, v, "%g" );
		FILE* foutS = fopen ( ".\\Results\\sMatrix.txt", "wb");
		gsl_vector_fprintf( foutS, s, "%g" );   
		cout<<"Done"<<endl;
	}
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
double pod::InterpolateCoefficients( const double& dx, const double& dy, 
      const double& dz, const int& row, gsl_matrix* coeffs, 
         const vector<int>& interp_columns )
{		

	// trilinear interpolation of the POD coefficients to get values at val 
   //vector

	double interpVal;

	double p000 = gsl_matrix_get( coeffs, row, interp_columns[0] );
	double p100 = gsl_matrix_get( coeffs, row, interp_columns[1] );
	double p010 = gsl_matrix_get( coeffs, row, interp_columns[2] );
	double p110 = gsl_matrix_get( coeffs, row, interp_columns[3] );
	double p001 = gsl_matrix_get( coeffs, row, interp_columns[4] );
	double p101 = gsl_matrix_get( coeffs, row, interp_columns[5] );
	double p011 = gsl_matrix_get( coeffs, row, interp_columns[6] );
	double p111 = gsl_matrix_get( coeffs, row, interp_columns[7] );

	double c0 = p000; 
	double c1 = (p100 -p000); 
	double c2 = (p010 -p000);
	double c3 = (p001 -p000); 
	double c4 = (p110 - p010 -p100 +p000);
	double c5 = (p011 -p001 - p010 +p000); 
	double c6 = (p101 -p001 -p100 +p000);
	double c7 = (p111 -p011 - p101 -p110 +p100 +p001 +p010 - p000);

	interpVal = c0 +c1*dx + c2*dy + c3*dz +c4*dx*dy + c5*dy*dz +c6*dz*dx + 
         c7*dx*dy*dz;

	return interpVal;
}

////////////////////////////////////////////////////////////////////////////////
void pod::search_range( const vector<double>& param, const double& val, 
         vector<int>& param_range )
{	

	// search monotonic data vector to find bounding values
	unsigned int i;
	for ( i=0;i<param.size()-1;i++ )
	{
		//search for the closest values
		if ( ( val >= param[i] ) && ( val <= param[i+1] ) )
			break;
	}
	param_range.push_back(i);
}




////////////////////////////////////
//////PUBLIC MEMBER FUNCTIONS///////
////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
void pod::ReadROM( )
{
	//USE THIS ONLY TO READ INFORMATION FROM THE INFO.TXT FILE
	//(not for input files)
	// read ROM information file
	ifstream fin;
	if ( DBG_POD )
		cout<<"Reading info.txt"<<endl;

	fin.open( ".\\Results\\info.txt" );
   //==============================================================================
	// read parameter vectors
   string firstLine;
   std::getline( fin, firstLine );
   //fin.getline(
   if ( DBG_POD > 2 )
      cout<<firstLine<<" Length is : "<<firstLine.length()<<endl;
   //grab number of parameters from this string
   
   string::iterator pos;
   int tmpCount = 0;
   for ( pos=firstLine.begin();pos!=firstLine.end();++pos )
   {
      if ( ( *pos ) == '\t' )
         tmpCount++;
   }
   //param_sizes.resize( tmpCount ); 
   param_sizes.resize( tmpCount + 1 ); //one parameter more than the number of spaces
   fin.seekg( 0, std::ios::beg );//reset file pointer

   //===========================================================================
   
   int tmpSize;   //temporary size var
   tmpSize = 0;
   for ( unsigned int i=0;i<param_sizes.size();i++ )
   {
      fin >> tmpSize;
      param_sizes.at(i) = tmpSize;
      if ( DBG_POD > 2 )
         cout<<"Param : "<<i<<"\t"<<tmpSize<<endl;
	  //fin >> tmpSize;
   }
   /*string line;
   for (unsigned int j=0;j<num_dims;j++)
   {*-
	   std::getline(fin, line);
   }*/
   //create a map to store the input parameter values
   double tmp;
   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
   {
	   vector< double > parameterHolder;
	   for ( int j=0;j<param_sizes.at(i);j++ )
	   {
		   fin >> tmp;
		   parameterHolder.push_back( tmp );
	   }
	   inputParameterMap.insert( make_pair( i, parameterHolder ) );
	   parameterHolder.clear( );
   }
   cout<<"INDEX================PARAMETERS====================="<<endl;
   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
   {

	   vector< double > parameterHolder;
	   parameterHolder = inputParameterMap.find(i)->second;
	   cout<<i<<"\t | \t";
	   for ( int j=0;j<param_sizes.at(i);j++ )
	   {
		   cout<<parameterHolder.at(j)<<"\t";
	   }
	   cout<<endl;
	   parameterHolder.clear( );
   }
   cout<<"===================================================="<<endl;
	// now read body of data
	int numEnergeticModes;
	fin >> numEnergeticModes;  // number of energetic modes
	fin >> sizeOfDataSet;      // read in size of data set, m, number of rows
	fin >> numDataSets;        // read in number of data sets,n, num of columns

	energeticModes.resize(numEnergeticModes);  // hack for now

   int totalNumDataSets = 1;   //var for assertion
   
   for ( unsigned int i=0;i<param_sizes.size();i++ )
   {
      totalNumDataSets = totalNumDataSets*param_sizes.at(i);
   }
   /*TODO: not sure what i was trying to do here
   if ( DBG_POD > 2 )
      cout<<"Total number of datasets : "<<totalNumDataSets<<endl;
	assert(numDataSets == totalNumDataSets );*/

	if ( DBG_POD ) cout<<"sizeOfDataSet :"<<sizeOfDataSet<<endl;
	if ( DBG_POD ) cout<<"numDataSets :"<<numDataSets<<endl;

	fin.close();

	// allocate matrices
	coeffs           = gsl_matrix_calloc( numDataSets, numDataSets );
	meanRemovedData  = gsl_matrix_calloc( sizeOfDataSet, numDataSets );
	vectorOfMeanData = gsl_vector_calloc( sizeOfDataSet );

	// read matrices	
	FILE *fptr;

	fptr=fopen(".\\Results\\coeffs.txt","rb");
	if(fptr==NULL)
	{
		cout<<"Unable to open coeffs.txt\n";
		exit(0);
	}

	if(gsl_matrix_fscanf (fptr, coeffs)==GSL_EFAILED)
	{
		cout<<"Error reading coeffs matrix\n";
		exit(0);
	}
	fclose(fptr);

	/////////////////////////////////////
	fptr=fopen(".\\Results\\uMatrix.txt","rb");
	if(fptr==NULL)
	{
		cout<<"Unable to open uMatrix.txt\n";
		exit(0);
	}

	if(gsl_matrix_fscanf (fptr, meanRemovedData)==GSL_EFAILED)
	{
		cout<<"Error reading u matrix\n";
		exit(0);
	}
	fclose(fptr);

	/////////////////////////////////////
	fptr=fopen(".\\Results\\mean.txt","rb");
	if(fptr==NULL)
	{
		cout<<"Unable to open mean.txt\n";
		exit(0);
	}

	if(gsl_vector_fscanf (fptr, vectorOfMeanData)==GSL_EFAILED)
	{
		cout<<"Error reading mean_data matrix\n";
		exit(0);
	}
	fclose(fptr);

}

////////////////////////////////////////////////////////////////////////////////
void pod::ReadDataFromFile( string fileName )
{
	//USE THIS ONLY TO READ ENSEMBLE DATA
	//Assuming a simple ASCII text file input format for the ensemble data. 
	//A translator will have to be written to get data from other file formats, 
   //Fluent/StarCD/VTK etc. 

	double tmp;    // temporary data var
	int counter = 0;
	int tmpSize;   //temporary size var

	ifstream fin;
	if ( DBG_POD )
		cout<<"Reading "<<fileName<<"...Done"<<endl;
	fin.open( fileName.c_str( ) );
   //==============================================================================
	// read parameter vectors
   string firstLine;
   getline( fin, firstLine );
   if ( DBG_POD > 2 )
      cout<<firstLine<<" Length is : "<<firstLine.length()<<endl;
   //grab number of parameters from this string
   
   string::iterator pos;
   int tmpCount = 0;
   for ( pos=firstLine.begin();pos!=firstLine.end();++pos )
   {
      if ( ( *pos ) == ' ' )
         tmpCount++;
   }
   
   param_sizes.resize( tmpCount + 1 ); //one parameter more than the number of spaces
   fin.seekg( 0, std::ios::beg );//reset file pointer
   
   //==============================================================================
   for ( unsigned int i=0;i<param_sizes.size();i++ )
   {
      fin >> tmpSize;
      param_sizes.at(i) = tmpSize;
      if ( DBG_POD > 2 )
         cout<<"Param : "<<i<<"\t"<<tmpSize<<endl;
   }
   
   //create a map to store the input parameter values
   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
   {
      vector< double > parameterHolder;
      for ( int j=0;j<param_sizes.at(i);j++ )
      {
         fin >> tmp;
         parameterHolder.push_back( tmp );
      }
      inputParameterMap.insert( make_pair( i, parameterHolder ) );
      parameterHolder.clear( );
   }
   cout<<"INDEX================PARAMETERS====================="<<endl;
   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
   {
      
      vector< double > parameterHolder;
      parameterHolder = inputParameterMap.find(i)->second;
      cout<<i<<"\t | \t";
      for ( int j=0;j<param_sizes.at(i);j++ )
      {
         cout<<parameterHolder.at(j)<<"\t";
      }
      cout<<endl;
      parameterHolder.clear( );
   }
   cout<<"===================================================="<<endl;

	// now read body of data
	//int numEnergeticModes;
	//fin >> numEnergeticModes;  // number of energetic modes
	fin >> sizeOfDataSet;   //read in size of data set, m, number of rows
	fin >> numDataSets;     //read in number of data sets,n, number of columns
   
   int totalNumDataSets = 1;   //var for assertion
   for ( unsigned int i=0;i<param_sizes.size();i++ )
   {
	  //TODO: not sure why this is here
      //totalNumDataSets = totalNumDataSets*param_sizes.at(i);
   }
   //if ( DBG_POD > 2 )
      //cout<<"Total number of datasets : "<<totalNumDataSets<<endl;
	//assert(numDataSets == totalNumDataSets );

	if ( DBG_POD ) cout<<"sizeOfDataSet :"<<sizeOfDataSet<<endl;
	if ( DBG_POD ) cout<<"numDataSets :"<<numDataSets<<endl;

	//allocate memory for matrix
	inputData = gsl_matrix_calloc( sizeOfDataSet, numDataSets );
	meanRemovedData = gsl_matrix_calloc( sizeOfDataSet, numDataSets );
	copyOfMeanRemovedData = gsl_matrix_calloc( sizeOfDataSet, numDataSets );

	//allocate memory for v matrix, nXn
	v = gsl_matrix_calloc( numDataSets, numDataSets );

	//allocate memory for s vector, nXn, diagonal vector
	s = gsl_vector_calloc( numDataSets );

	//read data into matrix
	while ( counter<sizeOfDataSet )
	{      
		//if ( DBG_POD ) cout<<"Counter is :"<<counter<<endl;      
		for ( int i=0;i<numDataSets;i++ )
		{
			fin >> tmp;         
			gsl_matrix_set( inputData, counter, i, tmp );
		}
		counter++;
	}   
	fin.close( );
}

////////////////////////////////////////////////////////////////////////////////
void pod::ComputeCoefficients( )
{
	//perform POD computations and return POD coefficients
	if ( DBG_POD )
		cout<<"Computing coefficients ..."<<endl;

	SubtractMeanOfData( );
	ComputeSVD(  );

	//ComputeSVD_GPU();
	//project left singular vectors onto original data to get POD coefficients

	if ( DBG_POD > 1 )
		cout<<"Assertion checked!!!"<<endl;

	//now project
	coeffs = gsl_matrix_calloc( numDataSets, numDataSets );
	if ( DBG_POD )
   {
      
      cout<<"projecting"<<"\t";
      cout<<"Size of coeffs : "<<coeffs->size1<<" "<<coeffs->size2<<endl;
   }

	if ( DBG_POD ) cout<<"Size of copyOfMeanRemovedData :"
      
      <<copyOfMeanRemovedData->size1<<" "<<copyOfMeanRemovedData->size2<<endl;
	if ( DBG_POD ) cout<<"Size of meanRemovedData :"<<meanRemovedData->size1
      <<" "<<meanRemovedData->size2<<endl;

	gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1.0, meanRemovedData, 
                  copyOfMeanRemovedData, 0.0, coeffs );
	FILE* foutCoeffs = fopen ( "./Results/coeffs.txt", "wb");
	gsl_matrix_fprintf( foutCoeffs, coeffs, "%g" );
}

////////////////////////////////////////////////////////////////////////////////
void pod::ComputeEnergy( double threshold )
{
	//compute the most energetic modes in the expansion, based on the threshold

	//all modes with energy levels <= threshold will not be considered 
   //in the POD expansion   

	double* energy;
	energy = new double [ (int)s->size ];
	//compute total energy
	double totalEnergy;
	totalEnergy = (double)(0);
	for ( unsigned int i=0;i<s->size;i++ )
	{
		totalEnergy = totalEnergy + gsl_vector_get( s, i );
	}
	if ( DBG_POD ) cout<<"Total Energy :"<<totalEnergy<<endl;

	for ( unsigned int i=0;i<s->size;i++ )
	{
		energy[i] = gsl_vector_get( s, i )/totalEnergy;
		if ( DBG_POD > 1 ) cout<<"Energy :"<<i<<" "<<energy[i]<<endl;
	}
   
   //now traverse the modes and find the most energetic modes greater than the 
   
   //threshold
	energeticModes.clear( );
	for ( unsigned int i=0;i<s->size;i++ )

	{
		if ( energy[i] >= threshold )
		{
			energeticModes.push_back( energy[i] );
		}
	}
	if ( DBG_POD )
	{
		cout.precision(5);
		cout<<(1.0-threshold)*(double)(100)<<
            " % Variation in data is captured by "<<
            energeticModes.size( )<<" modes"<<endl;
	}
	delete [] energy;
}
////////////////////////////////////////////////////////////////////////////////
void pod::SetInterpolationType( int type )
{
   interpolationType = type;
}
////////////////////////////////////////////////////////////////////////////////
void pod::ComputePODApproximation( const vector<double>& val, double, double )
{

	//--------->range checking for val to be implemented
	//to compute the POD approximation at the given val (the value of c)
	//cVals is the vector of values in which the parameter c varies
	//to interpolate coefficients loop over only the number of modes above the 
   //threshold

	interpolatedCoeffs.clear( );

	/*vector<int> param_range_low;
	param_range_low.clear();

	// find index range low for evaluation point
	search_range( param1, val[0], param_range_low );
	search_range( param2, val[1], param_range_low );
	search_range( param3, val[2], param_range_low );

	if ( DBG_POD > 1 ) 
	{
		cout<<"param1_range :"<<param_range_low[0]<<endl;
		cout<<"param2_range :"<<param_range_low[1]<<endl;
		cout<<"param3_range :"<<param_range_low[2]<<endl;
	}

	// find interpolation columns
	vector<int> interp_columns;

	for(int a=0;a<=1;a++)
	{
		for(int b=0;b<=1;b++)
		{
			for(int c=0;c<=1;c++)
			{
				interp_columns.push_back((param_range_low[0]+c)*
                  (param2.size()*param3.size())
					+(param_range_low[1]+b)*param3.size()+param_range_low[2]+a);
			}
		}
	}

	// find delta values for interpolation
	double dx = (val[0] - param1[param_range_low[0]])/
      (param1[param_range_low[0]+1] - param1[param_range_low[0]]); 
	double dy = (val[1] - param2[param_range_low[1]])/
      (param2[param_range_low[1]+1] - param2[param_range_low[1]]); 
	double dz = (val[2] - param3[param_range_low[2]])/
      (param3[param_range_low[2]+1] - param3[param_range_low[2]]); */
      
      
   //perform interpolation to get interpolated coefficients
   //TODO: not sure why i have this assert
	//assert( val.size() == inputParameterMap.size() );
   //set output multi_array, assuming a 4D array for now, 
   //will figure out n-D using boost later
   
   
   if ( DBG_POD )
      cout<<"Parameter sizes length : "<<param_sizes.size()<<endl;
   int counter = 0;
   
   double tmpCoeff;
   double tmp;
   if ( param_sizes.at(0) == 1 )
   {
	   vector<double> out;
	   vector<double> in;
	   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
	   {

		   vector< double > parameterHolder;
		   parameterHolder = inputParameterMap.find(i)->second;
		   cout<<i<<"\t | \t";
		   for ( int j=0;j<param_sizes.at(i);j++ )
		   {
			   tmp = parameterHolder.at(j);
			   in.push_back( tmp );
		   }
		   parameterHolder.clear( );
	   }
	   //in = inputParameterMap.f;
	   for ( unsigned int modes=0;modes<energeticModes.size();modes++ )
	   {
		   for ( int i=0;i<numDataSets;i++ )
		   {
			   out.push_back(gsl_matrix_get(coeffs, modes, i ));
		   }
		   	::Interpolation* inter = new ::Interpolation();
			inter->SetInputPoints1D( in, out );
			tmpCoeff = inter->EvaluateCubicSplineAt( val.at(0) );
			interpolatedCoeffs.push_back( tmpCoeff );
			out.clear();
			delete inter;
	   }
   }
   else if ( param_sizes.at(0) == 2 )
   {
	   vector<double> x;
	   vector<double> y;
	   gsl_matrix* z = gsl_matrix_calloc( energeticModes.size(), numDataSets );
	   for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
	   {
		   vector< double > parameterHolder;
		   parameterHolder = inputParameterMap.find(i)->second;
		   cout<<i<<"\t | \t";
		   tmp = parameterHolder.at(0);
		   x.push_back( tmp );
		   tmp = parameterHolder.at(1);
		   y.push_back( tmp );
		   parameterHolder.clear( );
	   }
	   for ( unsigned int modes=0;modes<energeticModes.size();modes++ )
	   {
		   for ( int i=0;i<numDataSets;i++ )
		   {
			   //out.push_back(gsl_matrix_get(coeffs, modes, i ));
			   gsl_matrix_set( z, modes, i, gsl_matrix_get(coeffs, modes, i ));
		   }
		   	/*::Interpolation* inter = new ::Interpolation();
			inter->SetInputPoints1D( in, out );
			tmpCoeff = inter->EvaluateCubicSplineAt( val.at(0) );
			interpolatedCoeffs.push_back( tmpCoeff );
			out.clear();
			delete inter;*/
		   ::BiCubicInterpolation* inter = new ::BiCubicInterpolation();
		   inter->SetInputPoints( x, y, coeffs );
		   //tmpCoeff = inter->EvaluateCubicSplineAt( val.at(0) );
		   tmpCoeff = inter->EvaluateAt( val.at(0), val.at(1) );
		   interpolatedCoeffs.push_back( tmpCoeff );
		   delete inter;
	   }
	   gsl_matrix_free( z );
   }
   else if ( param_sizes.at(0) == 3 )
   {
      typedef boost::multi_array<double, 3> array_type3;
      array_type3 outputMap3D(boost::extents
            
            [param_sizes.at(0)][param_sizes.at(1)][param_sizes.at(2)]);
      for ( unsigned int modes=0;modes<energeticModes.size();modes++ )
      {
         for ( int i=0;i<param_sizes.at(0);i++ )
         {
            
            for ( int j=0;j<param_sizes.at(1);j++ )
            {
               for ( int k=0;k<param_sizes.at(2);k++ )
               {
                  outputMap3D[i][j][k] = gsl_matrix_get( coeffs, modes, counter );
                  counter++;
               }
            }
         }
         
         tmpCoeff = this->Interpolation( val, outputMap3D );
         //cout<<"computed tmpCoeff : "<<tmpCoeff<<endl;
         interpolatedCoeffs.push_back( tmpCoeff );
         counter = 0;
      }  //modes loop
   }
   else if ( param_sizes.at(0) == 4 )
   {
      typedef boost::multi_array<double, 4> array_type4;
      array_type4 outputMap4D(boost::extents [param_sizes.at(0)]
            
            [param_sizes.at(1)][param_sizes.at(2)][param_sizes.at(3)]);
      
      counter = 0;
      for ( unsigned int modes=0;modes<energeticModes.size();modes++ )
      {
         for ( int i=0;i<param_sizes.at(0);i++ )
         {
            for ( int j=0;j<param_sizes.at(1);j++ )
            {
               for ( int k=0;k<param_sizes.at(2);k++ )
               {              
                  for ( int l=0;l<param_sizes.at(3);l++ )
                  {
                     outputMap4D[i][j][k][l] = gsl_matrix_get( coeffs, modes, 
                           
                           counter );
                     counter++;
                  }
               }
            }
         }
         tmpCoeff = this->Interpolation( val, outputMap4D );
         interpolatedCoeffs.push_back( tmpCoeff );
         counter = 0;
      }//modes loop
   }
   else if ( param_sizes.size() == 5 )
   {
	   typedef boost::multi_array<double, 5> array_type5;
	   array_type5 outputMap5D(boost::extents[param_sizes.at(0)]
	   [param_sizes.at(1)][param_sizes.at(2)][param_sizes.at(3)][param_sizes.at(4)]);
	   counter = 0;
      for ( unsigned int modes=0;modes<energeticModes.size();modes++ )
      {
         for ( int i=0;i<param_sizes.at(0);i++ )
         {
            for ( int j=0;j<param_sizes.at(1);j++ )
            {
               for ( int k=0;k<param_sizes.at(2);k++ )
               {
                  
                  for ( int l=0;l<param_sizes.at(3);l++ )
                  {
					  for ( int m=0;m<param_sizes.at(4);m++ )
					  {
						  outputMap5D[i][j][k][l][m] = gsl_matrix_get( coeffs, modes, counter );
						  counter++;
					  }
                  }
               }
            }
         }
         tmpCoeff = this->Interpolation( val, outputMap5D );
         interpolatedCoeffs.push_back( tmpCoeff );
         counter = 0;
      }//modes loop

   }
   if ( DBG_POD )
      cout<<" Finished interpolating all coefficients  : "<<
            interpolatedCoeffs.size()<<endl;

	//allocate memory
	PODapproximation = gsl_vector_calloc( sizeOfDataSet );

	//project the coefficients onto the left singular vectors
	//TODO : THIS CAN BE IMPROVED BY USING OpenMP
	time_t start, end;
	double tmpvar;
	time(&start);
	tmpvar = 0;
	for ( int j=0;j<sizeOfDataSet;j++ )
	{
		tmpvar = 0;
		for ( size_t i=0;i<energeticModes.size();i++ )
		{
         
         //cout<<"Interpolating : "<<i<<endl;
			tmpvar = tmpvar + interpolatedCoeffs.at(i)*
               gsl_matrix_get( meanRemovedData, j, i );
			gsl_vector_set( PODapproximation, j, tmpvar );
		}
	}

	//add the mean of the ensemble back to the projection   
	ofstream foutROM;
	foutROM.open( ".\\Results\\ROM.txt" );
	for ( int j=0;j<sizeOfDataSet;j++ )
	{
		tmpvar = gsl_vector_get( PODapproximation, j );
		tmpvar = tmpvar + gsl_vector_get( vectorOfMeanData, j );
		gsl_vector_set( PODapproximation, j, tmpvar );

		foutROM<<gsl_vector_get( PODapproximation, j )<<endl;      
		tmpvar = 0;
	}
	foutROM.close( );
	time(&end);
	cout<<"Time for POD computation :"<<difftime(end, start)<<endl;
}

////////////////////////////////////////////////////////////////////////////////
void pod::ComputePODApproximation( vector < double> cVals, int numModes )
{
	//--------->range checking for val to be implemented
	//to compute the POD approximation at the given val (the value of c)
	//cVals is the vector of values in which the parameter c varies
	//to interpolate coefficients loop over only the number of modes above the threshold
	int cl=-1, cu = -1;
	vector<double> thisvec;
	vector<double> nextvec;
	double r, rl, ru;

	for ( size_t i =0;i< inputParameterMap.size()-1 ; i++ )
	{
		thisvec = inputParameterMap.find( i )->second;
		nextvec = inputParameterMap.find( i+1 )->second;

		for ( size_t i=0;i<thisvec.size();i++ )
			cout<<thisvec.at( i )<<" ";
		if ( ( cVals > thisvec ) && ( cVals < nextvec ) )
		{
			cl = i;
			cu = i+1;
			cout<<"found interp points "<<endl;
			//cout<<thisvec<<"\t"<<nextvec<<endl;
			break;
		}
		else if ( ( cVals == thisvec ) || ( cVals == nextvec ) )
		{
			cout<<"found point in phase space"<<endl;
			break;
			//assign cl and cu
		}

		cout<<endl;
	}
	if ( (cl<0) && (cu<0) )
	{
		cout<<"Input vector not in phase space"<<endl;
		exit( -1 );
	}
	interpolatedCoeffs.clear( );
	double tmp; //temporary data var

	unsigned int i;
	r = length( cVals );
	rl = length( thisvec );
	ru = length( nextvec );
	//SOME CODE TO HELP WITH DEBUGGING, IF NEEDED
	/*for ( i=0;i<cVals.size()-1;i++ )
	{
		//search for the closest values
		if ( ( val >= cVals.at(i) ) && ( val <= cVals.at(i+1) ) )
			break;
	}

	cl = i;   
	cu = i+1;
	if ( DBG_POD > 1 ) cout<<"Range found at :"<<cl<<" "<<cu<<endl;*/
	if ( numModes == 0 )
		numModes = energeticModes.size();

	for ( size_t i=0;i<numModes;i++ )
	{
		//SOME CODE TO HELP WITH DEBUGGING, IF NEEDED
		/*tmp = gsl_matrix_get( coeffs, i, cl ) + 
			( gsl_matrix_get( coeffs, i, cu ) - gsl_matrix_get( coeffs, i, cl ) )*
			( val - cVals.at(cl) )/( cVals.at(cu) - cVals.at(cl) );*/
		tmp = gsl_matrix_get( coeffs, i, cl ) + 
			( gsl_matrix_get( coeffs, i, cu ) - gsl_matrix_get( coeffs, i, cl ) )*
			( r - rl )/( ru - rl );
		interpolatedCoeffs.push_back( tmp );
		if ( DBG_POD ) cout<<"Interpolated coefficients :"<<interpolatedCoeffs.at(i)<<endl;
	}

	//allocate memory
	PODapproximation = gsl_vector_calloc( sizeOfDataSet );
	//project the coefficients onto the left singular vectors

	double tmpvar;
	tmpvar = 0;
	for ( int j=0;j<sizeOfDataSet;j++ )
	{
		tmpvar = 0;

		for ( size_t i=0;i<numModes;i++ )
		{
			tmpvar = tmpvar + interpolatedCoeffs.at(i)*gsl_matrix_get( meanRemovedData, j, i );
			gsl_vector_set( PODapproximation, j, tmpvar );
		}
	}

	//add the mean of the ensemble back to the projection

	ofstream foutROM;
	foutROM.open( "./Results/ROM.txt" );
	for ( int j=0;j<sizeOfDataSet;j++ )
	{
		tmpvar = gsl_vector_get( PODapproximation, j );
		tmpvar = tmpvar + gsl_vector_get( vectorOfMeanData, j );
		gsl_vector_set( PODapproximation, j, tmpvar );

		foutROM<<gsl_vector_get( PODapproximation, j )<<endl;
		tmpvar = 0;
	}
	foutROM.close( );
}

////////////////////////////////////////////////////////////////////////////////
gsl_vector* pod::GetEigenVector( int index )
{
	//get the eigenvector corresponding to the user input index 
	//can be used to plot and analyse the flow patterns/data
	//index is the column number
	double tmp;
	eigenVector = gsl_vector_calloc( numDataSets );
	for ( int i=0;i<numDataSets;i++ )
	{
		tmp = gsl_matrix_get( v, i, index );
		gsl_vector_set( eigenVector, i, tmp );
	}
	if ( DBG_POD == 1 )
		cout<<"Returning eigenvector :"<<index<<endl;

	return eigenVector;
}

////////////////////////////////////////////////////////////////////////////////
int pod::GetEnergeticModes( )
{
	return energeticModes.size();
}
void pod::WriteMode( int index )
{
	//get the ith modes, i = index
	//the modes can be used to analyze the flow/data patterns

	//gives an idea for which mode contributes to the kind of pattern in the data
	//writes the mode to a file in the Results directory
	string fName = "./Results/Mode";
	ostringstream fn;
	fn<<fName<<index<<".txt";
	//allocate memory for mode   
	gsl_vector* PODMode = gsl_vector_calloc( sizeOfDataSet );
	double tmpvar;


	tmpvar = 0;
	if ( index >= (int) interpolatedCoeffs.size() )
		cout<<"Mode "<<index<<" is not within the energy limits set, will not"
		<< " compute this mode"<<endl;
	if ( index < (int) interpolatedCoeffs.size() )
	{
		cout<<"Writing..."<<fn.str().c_str()<<endl;
		for ( int j=0;j<sizeOfDataSet;j++ )
		{
			tmpvar = interpolatedCoeffs.at(index)*
				gsl_matrix_get( meanRemovedData, j, index );
			gsl_vector_set( PODMode, j, tmpvar );
		}

		FILE* foutMode = fopen ( fn.str().c_str(), "wb");
		gsl_vector_fprintf( foutMode, PODMode, "%g" );
	}
	gsl_vector_free( PODMode );

}

////////////////////////////////////////////////////////////////////////////////
void pod::WriteInfo( )
{

	// read ROM information file
	ofstream fout;
	if ( DBG_POD )
		cout<<"Writing info.txt"<<endl;
	fout.open( ".\\Results\\info.txt" );

	// write parameter vectors

	/*fout << param1.size()<<" "<<param2.size()<<" "<<param3.size()<<endl;

	for(unsigned int i=0;i<param1.size();i++)
	{
		fout<<param1[i]<<" ";
	}
	fout<<endl;

	for(unsigned int i=0;i<param2.size();i++)
	{
		fout<<param2[i]<<" ";
	}
	fout<<endl;

	for(unsigned int i=0;i<param3.size();i++)
	{
		fout<<param3[i]<<" ";
	}
	fout<<endl;*/
	for ( int i=0;i<param_sizes.size();i++)
	{
		fout << param_sizes.at(i);
		if (i!=param_sizes.size()-1)
			fout << "\t";
	}
	fout << endl;
	for ( unsigned int i=0;i<param_sizes.size();i++ )  //loop over num of parameters
	{

		vector< double > parameterHolder;
		parameterHolder = inputParameterMap.find(i)->second;
		for ( int j=0;j<param_sizes.at(i);j++ )
		{
			fout<<parameterHolder.at(j)<<"\t";
		}
		fout<<endl;
		parameterHolder.clear( );
	}

	// now write sizes
	fout << energeticModes.size()<<endl;  // number of energetic modes
	fout << sizeOfDataSet<<endl;  // read in size of data set, m, number of rows
	fout << numDataSets<<endl;     // read in num of data sets,n, num of columns

	fout.close();

}

////////////////////////////////////////////////////////////////////////////////
void pod::SetMode( const string& Mode)
{
	mode = Mode;
}
////////////////////////////////////////////////////////////////////////////////
